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ABSTRACT 

We present results from the first generation of multi-dimensional hydrodynamic 
core-collapse simulations in full general relativity (GR) that include an approximate 
treatment of neutrino transport. Using a Ml closure scheme with an analytic variable 
Eddington factor, we solve the energy-independent set of radiation energy and mo- 
mentum based on the Thome's momentum formalism. To simplify the source terms 
of the transport equations, a methodology of multiflavour neutrino leakage scheme is 
partly employed. Our newly developed code is designed to evolve the Einstein field 
equation together with the GR radiation hydrodynamic equations in a self-consistent 
manner while satisfying the Hamiltonian and momentum constraints. An adaptive- 
mesh-refinement technique implemented in the three-dimensional (3D) code enables us 
to follow the dynamics starting from the onset of gravitational core-collapse of a 15 Mq 
star, through bounce, up to about 100 ms postbounce in this study. By computing 
four models that differ according to ID to 3D and by switching from special relativistic 
(SR) to GR hydrodynamics, we study how the spacial multi-dimensionality and GR 
would affect the dynamics in the early postbounce phase. Our 3D results support the 
anticipation in previous ID results that the neutrino luminosity and average neutrino 
energy of any neutrino flavor in the postbounce phase increase when switching from SR 
to GR hydrodynamics. This is because the deeper gravitational well of GR produces 
more compact core structures, and thus hotter neutrino spheres at smaller radii. By 
analyzing the residency timescale to the neutrino-heating timescale in the gain region, 
we show that the criterion to initiate neutrino-driven explosions can be most easily sat- 
isfied in 3D models, irrespective of SR or GR hydrodynamics. Our results suggest that 
the combination of GR and 3D hydrodynamics provides the most favorable condition 
to drive a robust neutrino-driven explosion. To draw a solid conclusion, the energy and 
angle dependence of neutrino transport should be accurately incorporated in our full 
GR code with the use of more detailed set of weak interactions. This work is the very 
first step before going up the long and winding road. 

Subject headings: supernovae: collapse — neutrinos — hydrodynamics — general rela- 
tivity 
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1. Introduction 

Core-collapse supernova simulations have been counted as one of the most challenging subjects 
in computational astrophysics. The four fundamental forces of nature are all at play; the collapsing 
iron core bounces due to strong interactions; weak interactions determine the energy and lepton 
number loss in the core via the transport of neutrinos; electromagnetic interactions determine the 
properties of the stellar gas; general relativity plays an important role due to the compactness of the 
proto-neutron star and also due to high velocities of the collapsing material outside. Naturally, such 
physical richness ranging from a microphysical scale (i.e. femto- meter scale) of strong/weak inter- 
actions to a macrophysical scale of stellar explosions has long attracted the interest of researchers, 
necessitating a world-wide, multi-disciplinary collaboration to clarify the theory of massive stellar 
core-collapse and the formation mechanisms of compact objects. 

Ever since the first numerical simulation of such events (Colgate & White 1966), the neutrino- 
heating mechanism (Wilson 1985; Bethe &; Wilson 1985), in which a stalled bounce shock could 
be revived via neutrino absorption on a timescale of several hundred milliseconds after bounce, 
has been the working hypothesis of supernova theorists for these ~ 45 years. However, the sim- 
plest, spherically-symmetric (ID) form of this mechanism fails to blow up canonical massive stars 
(Rampp & Janka 2000; Liebendorfer et al. 2001; Thompson et al. 2003; Sumiyoshi et al. 2005). 
Pushed by mounting supernova observations of the blast morphology (e.g., Wang et al. (2001); 
Maeda et al. (2008); Tanaka et al. (2009), and references therein), it is now almost certain that 
the breaking of the spherical symmetry holds the key to solve the supernova problem. So far a 
number of multidimensional (multi-D) hydrodynamic simulations have shown that hydrodynamic 
motions associated with convectivc overturn (e.g., Herant et al. (1992); Burrows et al. (1995); Janka 
& Miiller (1996); Fryer et al. (2002); Fryer (2004)) and the Standing-Accretion-Shock-Instability 
(SASI, e.g., Blondin et al. (2003); Scheck et al. (2004); Scheck et al. (2006); Ohnishi et al. (2006); 
Ohnishi et al. (2007); Foghzzo et al. (2006); Iwakami et al. (2008); Iwakami et al. (2009); Murphy 
& Burrows (2008); Fernandez & Thompson (2009b, a), and references therein) can help the onset 
of the neutrino-driven explosion. 

In fact, the neutrino-driven explosions have been obtained in the following first-principle two- 
(2D) and three-(3D) dimensional simulations in which the spectral neutrino transport is solved at 
various levels of approximations (e.g., Miiller et al. (2011); Ott et al. (2011); Kotake (2011) for recent 
status reports). Among them are the 2D neutrino-radiation- hydrodynamic simulations by Buras 
et al. (2006a,b); Marek & Janka (2009) who included one of the best available neutrino transfer 
approximations by the ray-by-ray variable Eddington factor method, by Brucnn et al. (2010) who 
included a ray-by-ray multi-group flux-limited diffusion (MGFLD) transport with the best available 
weak interactions, and by Suwa et al. (2010, 2011) who employed a ray-by-ray isotropic diffusion 
source approximation (IDSA) (Liebendorfer et al. 2009) with a reduced set of weak interactions. By 
extending the 2D modules in Suwa et al. (2010), Takiwaki et al. (2011) recently reported neutrino- 
driven explosion models in 3D for an 11.2 Mq star. They pointed out whether 3D effects would 
help explosions or not is sensitive to the employed numerical resolutions (see also Hanke et al. 
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(2011); Nordhaus et al. (2010)). They argued that future peta- and exa-scale resources are at least 
needed to draw a robust conclusion of the 3D effects. 

In addition to the 3D effects, impacts of general relativity (GR) on the neutrino-driven mecha- 
nism stand out among the biggest open questions in the supernova theory. It should be remembered 
that using newly derived Einstein equations (Misner & Sharp 1964), the consideration of GR was 
standard in the pioneering era of supernova simulations (e.g.. May &; White (1966)). One year after 
Golgate & White (1966), Schwartz (1967) reported the first fully GR simulation of stellar collapse 
to study the supernova mechanism, who implemented a gray transport of neutrino diffusion in 
the ID GR hydrodynamics^. Using GR Boltzmann equations derived by Lindquist (1966), Wilson 
(1971) developed a ID GR-radiation-hydrodynamic code including a more realistic (at the time) 
description of the collisional term than the one in Schwartz (1967). By performing ID GR hydrody- 
namic simulations that included a leakage scheme for neutrino cooling, hydrodynamical properties 
up to the prompt shock stagnation were studied in detail (van Riper 1979; van Riper k, Lattimer 
1981; van Riper 1982). These pioneering studies, albeit using a much simplified neutrino physics 
than today, did provide a bottom-line of our current understanding of the supernova mechanism 
(see Bruenn et al. (2001) for a complete list of references for the early GR studies). In the middle 
of the 1980s, Bruenn (1985) developed a code that coupled ID GR hydrodynamics to the MGFLD 
transport of order (f/c) including the so-called standard set of neutrino interactions. Since the late 
1990s, the ultimate ID simulations, in which the GR Boltzmann transport is coupled to ID GR 
hydrodynamics, have been made feasible by Sumiyoshi-Yamada et al. (Yamada 1997; Yamada et al. 
1999; Sumiyoshi et al. 2005, 2007)^ and by Licbcndorfcr-Mczzacappa-Bruenn et al. (Mczzacappa 
& Matzner 1989; Bruenn et al. 2001; Liebendorfer et al. 2001; Liebendorfer et al. 2001, 2004) (and 
by their collaborators). 

Among them, Bruenn et al. (2001) presented evidence that average neutrino energy of any 
neutrino flavor during the shock reheating phase increase when switching from Newtonian to GR 
hydrodynamics. They also pointed out that the increase is larger in magnitude compared to the 
decrease due to redshift effects and gravitational time dilation. By employing the currently best 
available weak interactions, Lentz et al. (2011) very recently reported the update of Bruenn et al. 
(2001). They showed that the omission of observer corrections in the transport equation particularly 
does harm to drive the neutrino-driven explosions. In these full-fledged ID simulations, a commonly 
observed disadvantageous aspect of GR to drive neutrino-driven explosions is that the residency 
time of material in the gain region becomes shorter due to the stronger gravitational pull. As 
a result of these competing ingredients in the end, GR works disadvantageously to facilitate the 



^Citing from his paper, "In this calculation, the neutrino luminosity of the core is found to he IQ^^ erg/s, or 1/2 
a solar rest mass per second !! .... This is the mechanism which the supernova explodes". The neutrino luminosity 
rarely becomes so high in the modern simulations, but it is surprising that the potential impact of GR on the 
neutrino-heating mechanism was already indicated in the very first GR simulation. 

^Very recently, they reported their success to develop the first multi-angle, multi-energy neutrino transport code 
in 3D (Sumiyoshi & Yamada 2012). 
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neutrino-driven explosions in ID. In fact, the maximum shock extent in the postbounce phase is 
shown to be 20% smaller when switching from Newtonian to GR hydrodynamics (e.g., Figure 2 in 
Lentz et al. (2011)). 

Among the most up-to-date multi-D models with spectral neutrino transport mentioned earlier, 
GR effects are at best attempted to be modeled by using a modified gravitational potential that 
takes into account a ID, post-Newtonian correction (Buras et al. 2006a,b; Marek & Janka 2009; 
Bruenn et al. 2010). A possible drawback of this prescription is that a conservation law for the 
total energy cannot be guaranteed by adding an artificial term to the Poisson equation of self- 
gravity. Since the energy reservoir of the supernova engines is the gravitational binding energy, any 
potential inaccuracies in the argument of gravity would be better eliminated. There are a number 
of relativistic simulations of massive stellar collapse in full GR (e.g., 2D (Shibata &: Sekiguchi 
2005a) or 3D (Shibata & Sekiguchi 2005b; Ott et al. 2007), and references therein) or using the 
conformally-flatness approximation (CFC) (e.g., Dimmelmeier et al. (2002); Cordero-Carrion et 
al. (2009)). Although extensive attempts have been made to include microphysics such as by the 
Ye formula (Liebendorfer 2005) or by a neutrino leakage scheme (Sekiguchi 2010), the effects of 
neutrino heating have yet to be included in them, which is a main hindrance to study the GR 
effects on the multi-D neutrino-driven mechanism^. 

In this paper, we present a new fully GR code for multi-D hydrodynamic supernova simulations 
in which an approximate neutrino transport is implemented. The code is a marriage of an adaptive- 
mesh-refinement (AMR), conservative 3D GR magnetohydrodynamic (MHD) code developed by 
Kuroda Sz Umeda (2010), and the approximate neutrino transport code that we newly develop 
in this work. The spacetime treatment in our full GR code is based on the Baumgarte-Shapiro- 
Shibata-Nakamura (BSSN) formalism (see, e.g., Shibata &; Nakamura 1995; Baumgarte & Shapiro 
1999). The hydrodynamics can be solved either in full GR or in special relativity (SR), which allows 
us to investigate the GR effects on the supernova dynamics. Using a Ml closure scheme with an 
analytic variable Eddington factor, we solve the energy-independent set of radiation energy and 
momentum. This part is based on the partial implementation of the Thome's momentum formalism, 
which is recently extended by Shibata et al. (2011) in a more suitable manner applicable to the 
neutrino transport problem. Similar to the IDSA scheme (Liebendorfer et al. 2009) , we conceptually 
divide the neutrinos into two parts, which are "trapped" and "free-streaming" neutrinos. By doing 
so, we model the source terms of the transport equations to be expressed in a simplified manner 
with the use of a multi-flavor neutrino leakage scheme (e.g., Rosswog & Liebendorfer (2003)). 
Our newly developed code is designed to evolve the Einstein field equation together with the GR 
radiation hydrodynamic equations in a self-consistent manner while satisfying the Hamiltonian 
and momentum constraints. An adaptive-mesh-refinement technique implemented in the 3D code 
enables us to follow the dynamics starting from the onset of gravitational core-collapse of a 15 



^Very recently, Miiller et al. (2012) reported explosions for 11.2 and 15Mq stars based on their 2D GR simulations 
in CFC with detailed neutrino transport (Miiller et al. 2010) similar to Buras et al. (2006a). 
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Mq star, through bounce, up to about 100 ms postbounce in this study. For the 15 Mq star, the 
neutrino-driven explosions are expected to take place later than ~ 200 ms postbounce at the earliest 
(e.g., Bruenn et al. (2010); Marek & Janka (2009)). However it is computationally too expensive to 
follow such a long-term evolution in our full 3D GR simulations. Albeit limited to the rather early 
postbounce phase, we would exploratory study possible GR effects in the multi-D neutrino-driven 
mechanism by comparing ID to 3D results and by switching from SR to GR hydrodynamics. 

This paper is organized as follows: In section 2, after we introduce the model concept of 
the approximate GR transport scheme, we summarize the governing equations of hydrodynamics 
and neutrino transport in detail. The main results are presented in Section 3. We summarize 
our results and discuss their implications in Section 4. Note that geometrized units are used 
throughout Sections 2 to 3, i.e. both the speed of light and the gravitational constant are set to 
unity: G = c = 1. Greek indices run from to 3, Latin indices from 1 to 3. 



2. Basic Equations for General Relativistic Neutrino-Radiation Hydrodynamics 

Our newly developed code consists of the three parts, in which the evolution equations of 
metric, hydrodynamics, and neutrino radiation are solved, respectively. As will be mentioned, each 
of them is solved in an operator-splitting manner, but the system evolves self-consistently as a 
whole satisfying the Hamiltonian and momentum constraints. Before going into details, we shortly 
describe the bottom-line how to add radiation to GR hydrodynamics. 

The starting-point is the conservation of the total energy (fluid -|- radiation) , 

^ a-' (total) - ^ a-' (fluid) + ^""'(i/) " ^' V^-J 

where T'^^otai)' -^(fluid)' ^^'^ '^{f) '^^ stress-energy tensor of the total energy, fluid, and neutrino 
radiation, respectively. Then Equation (1) can be decomposed as, 

Var(td) = -Q'^ (2) 

and 

VaTj^ = (3) 

where represents the source terms that describe the exchange of energy and momentum between 
fluid and radiation. If would be given, it is rather straightforward to evolve Equation (2) 
following standard procedures in numerical relativity. Accordingly, what we focus on in this section 
is how to determine the source terms and the evolution equations of neutrinos (i.e. the concrete 
form of the left-hand-side of Equation (3)). In doing so, there will appear many terms related to 
GR such as e^^, di^K. etc. So, after we briefly summarize the BSSN formalism in the next section, 
we first present the transport equations in section 2.2 and then the source terms in section 2.3. 
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2.1. Metric equations 

We write the spacetime metric in the standard (3+1) form: 

ds'^ = -a^dt^ + -iij (dx' + P'dt) {dx^ + /3^' dt) , (4) 

where a, /3% and 7jj are the lapse, shift, and spatial metric, respectively. The extrinsic curvature 
Kij is defined by 

{dt - Ci3)-fij = -2aKij, (5) 
where Cj} is the Lie derivative with respect to The evolution of jij and Kij is governed by the 
Einstein equation Gfj^i, = SttT^,^ (tg^ai) ) where G^u is the Einstein tensor and (total) total 
stress-energy tensor (e.g.. Equation (1)). 

We evolve 7ij and Kij using the BSSN formulation (Baumgarte &: Shapiro 1999; Shibata & 
Sekiguchi 2005b; Duez et al. 2006), in which the fundamental variables are 

^ = ^ln[det(7i,)] , (6) 

lij ^ (7) 
K = j'^Kij, (8) 

Aij = e-''^{Kij-^jijK) , (9) 

r = (10) 

The Einstein equation gives rise to the evolution equations for the BSSN variables as, 

{dt-£p)jij = -2aAij (11) 
{dt-C0)4> = ~aK (12) 



trf 



{dt - Cp)Aij = e-""^ [a{Rij - 87r7i^7,vrJ^tai) " DiDja\ + a{KAij - 2AikTAji) 

(13) 

{dt-Cp)K = -Aa + a{A,A^^+Ky3)+4na{n^n,T^;i^^^+j'^ji^jj,T^^^^^ (14) 
{dt-P''dk)f = 167rf^7«,.n.Tj^tai) 

-2a{^rKj-6A'^<l,,j-f%A^'^) 

+f'p'jk + - ^'Pij + li^'p'j + p'i^ij - 2i'^«,„ (15) 

where D denotes covariant derivative operator associated with ^ij, A = L)*L'i,"trf" denotes the 
trace-free operator, = {—a, 0) is the time-like unit vector normal to the t = constant time slices. 
In Equation (13), the explicit form of DiDja reads 

DiDja = didja — 



didja 



f^j + 2 (^6pi<j> + 6fdj(l) - jijj'^'dicl)) ] dka. (16) 
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FoUowing Alcubierre & Briigmann (2001), the gauge is specified by the 1+log lapse, 

dta = P'dia - 2aK, (17) 

and by the Gamma-driver-shift, 

dtp' = kdtr\ (18) 

here we chose k = 1. For further information with code verification of the metric solver, see Kuroda 
&: Umeda (2010)^. Having summarized the bottom-line of the BSSN formalism, we are now ready 
to discuss the transport equations. 



2.2. Neutrino Transport Equations 

To determine the transport equations in GR, we follow the truncated momentum formalism 
(Thorne 1981), which is recently extended by Shibata et al. (2011) in a suitable form for the neutrino 
transport problem. The starting point is to define the radiation stress-energy tensor as, 

T^^ff ^ E^,)n^n^ + F^^fn^ + F^^fn^ + P^^f^, (19) 

where E(^^^, and P(,/): is the radiation energy, flux, pressure measured in the laboratory frame, 
respectively. Conversely, , F^^^ , and are given by the stress-energy tensor as. 

In the following, radiation variables are all defined in the laboratory frame unless otherwise stated. 

According to Shibata et al. (2011), the evolution equations of radiation energy and radiation 
flux in Equation (19) can be written as 

5t(e6'^^(,)) + diie'^'f'iaFi^^ - P'E^,))] = e'^'l' (aP'^ K^j - F^^^^a - aQ^n^), (21) 

and 

9t(e6'^F(,).)+a,[e6-^(aP(,)^-/3^F(,).)] = e6'^[-S(,)5,a+F(,) .S^^S^ + (a/2)Pg5,7,■fe+aQ^7^^], (22) 

where denotes the source terms. For the three radiation variables (-£'(;/), -^(V)' ■^(^)) Equations 
(21,22), we have only two sets of the equation. Here we employ the so-called Ml closure (Levermore 
1984), in which the radiation pressure is related to the radiation energy and flux by an analytical 
closure relation (i.e. P(,^)(£^(^), F(y))) as, 

p ij — "^X ~ ^ pij I 3(1 — x) jj 

P(u) - ^ Pthin + 2 ^ ^ 



*ln Kuroda & Umeda (2010), an auxiliary variable Fi = S^'^dkjij was evolved instead of f 
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where % represents the variable Eddington factor, P^^^^ and -P^'^ick corresponds to the radiation 
pressure in the optically thin and thick limit, respectively. 



For the variable Eddington factor %, we employ the one proposed by Levermore (1984), 

3 + 4^2 
5 + 2V4 - 3^2 ' 



X = ^ ^ (24) 



- (25) 

It can be readily checked that in the optically thick limit, x ~^ 1/3 because — >■ 0, then P*^ — >■ 
-^thick' while in the optically thin limit, x ^ 1 because 1, then P*^ -^thin- 

Following Audit et al. (2002); Shibata et al. (2011), the following forms of pH^^ and pH^^^ are 
adopted, 

pi pj 

^hin = Ej^, (26) 

and 

PiU = J ^^'+^^^V'W ^ yfe^.^^ ^ ^^k^i^^^ ^27) 

respectively. By this choice, the radiation flux naturally changes with radius (r) as ~ 1 /r^ in the 
low opacity regime (e.g., Appendix C). This may sound quite straightforward, but it is one of 
the most important issue for the purpose of this work, because the radiation neutrino flux in the 
semi-transparent regions holds the key to the success or failure of the neutrino-heating mechanism. 
J, T-L in Equation (27) denotes the Eddington moments in the comoving frame, which are related 
to those in the laboratory frame as, 

J = E(^,)W^ - 2WF^,)'ui + P^^fu^uj, (28) 

and 

= {E^,)W - F^,)'ui){n" - W^n") + l^/i^F(,)^ - KujP^.^'^, (29) 
where W = au^ is the Lorentz factor, is the projection operator deflned by 

haP = 3ap+UaUp. (30) 

Having summarized the closed set of the two-moment transport equations, we are now going to 
discuss the source terms (;,Q^) in the next section. 



2.3. Source terms 



To model the source terms, we follow the idea of the IDSA scheme (Liebendorfer et al. 2009), in 
which neutrinos are conceptually divided into two parts, which are "trapped" and "free-streaming" 
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neutrinos, respectively. We also utilize a methodology of multi-flavor neutrino leakage scheme (e.g., 
Rosswog & Liebendorfer (2003)) to simplify the description of the source terms. 

Figure 1 illustrates the concept how to estimate the source terms. To describe the neutrino- 
matter coupling, we need to ask at least three actors, namely "matter", "trapped neutrino", and 
"streaming neutrino" , to appear in the playground of the supernova core. The neutrino sphere is an 
important quantity to describe the relationship between them. The position of the neutrino sphere^ 
at which the neutrino optical depth comes close to unity, is indicated by r,^ = 2/3 in Figure 1. The 
trapped neutrinos are always coupled with matter (through /3-equilibrium) and their temperature 
is the same as that of matter. On the other hand, temperature of the free-streaming neutrinos 
cannot be determined locally owing to its decoupling nature from matter, which was the reason 
that we have to solve the evolution equations. So these two represent a two extreme limit. 

In Figure 1, trapped neutrinos are denoted by "t'trap" (in a diamond shape colored by grey), 
which is illustrated to shake hands with matter inside the neutrino sphere (inside the region enclosed 
by Tjy = 2/3). There the trapped neutrinos dominate over the streaming neutrinos (denoted by 
"i^stream" (in a jaggy circle colored by orange) in Figure 1), which is vice versa outside the neutrino 
sphere. This is illustrated in such a way that "i^trap" is bigger than "I'stream" inside the neutrino 
sphere, which is vice versa outside the neutrino sphere. For "I'stream" outside the neutrino sphere, 
the jaggy circle is drawn to have several tails, by which we intend to express that it can travel much 
more freely in the free-streaming regime. 

In Figure 1, the three actors are connected by thick arrows (in blue or red), each of them is 
labeled by Q^^j (in blue), Qf^^ (in blue), or Q^'^ (in red), representing the couplings in-between. 
The arrows colored by blue {Q^jjp Q^r) represent neutrino cooling, while the arrow in red {Q^'^) 
does neutrino heating. The neutrino cooling means that energy is transferred from matter (or from 
trapped neutrinos) to streaming neutrinos that carry the imparted energy away from the system. 
On the other hand, the neutrino heating proceeds by energy transfer from streaming neutrinos to 
matter (see Qi^'^ in Figure 1). Finally Q^l^f in Figure 1 represents the cooling by neutrinos leaking 
out from opaque regions inside the neutrino sphere by diffusion.^ 

Looking at Figure 1 again, the source term of Equations (21,22) can be readily defined as. 



where each of the cooling {Q^' ) and heating(Q'^' ) term is calculated in the present scheme as. 



®here defined for the average neutrino energy for simplicity, 

®Note inside the neutrino sphere, neutrino heating locally balances with neutrino cooling by weak interactions due 
to /3-equilibrium. So as a net, the diffusion-mediated cooling becomes dominant there. 



(31) 





(32) 



-10- 



Q'"'" = J2 (^~^"^"elk,{-Jv'^''-'Hi)- (33) 

Before we go into details, we need to draw a caution that wc introduced the concept of the streaming 
and trapped neutrinos only for the sake of (better) explanation of our approximate treatment. 
Actually the sum of the trapped and streaming is transported by Equations (21,22) with the source 
terms described above. We design the source terms in such a way to connect the heating/cooling 
terms smoothly between the diffusion and free-streaming limit, which is basically similar to the 
concept of the Ml closure relation. 

The cooling term (Q^'*^) consists of Qu,diff and Q^[intr^ which accounts for neutrino coohng 
by diffusion out of the neutrino sphere and the one determined locally outside the neutrino sphere, 
respectively (e.g.. Figure 1). Following van Riper & Lattimer (1981), the terms of 1 — e~^'''^'' and 
^-l3vT„ appearing in Equation (32) are introduced to smoothly connect the two quantities {Qu^diff 
and Q^yintr) foi" the semi-transparent regime. Here Ty represents the optical depth of neutrinos and 
I3i, is a model parameter that we determine by the comparison with a spectral neutrino transport 
calculation (see Appendix A.l). With these terms bridging the two limits, it is easy to see that Q'*'^ 
approaches to Qu,diff for the diffusion limit {Ty — > oo), and it does to Q^^intr free-streaming 
limit (tj, -)• 0). 

Following the neutrino leakage scheme (e.g., Epstein & Pethick (1981); van Riper & Lattimer 
(1981); Kotake et al. (2003)), the diffusion cooling rate (Qy^diff) can be given as 

/£ 71/ i S ) 
LJf, \ dsi., (34) 

where the right-hand-side of Equation (34) is simply proportional to the leakage of the neutrino 
energy density : Uy [erg/cm^] divided by the diffusion timescale Ty'^^^ [s]. More details to estimate 
these quantities as well as how to determine the neutrino sphere are summarized in Appendix A. 

Charged Current Interactions 

pi>e ^ e^n 
VeA O e~A' 
Neutral Current Interactions 
up -f-> vp 

vn -H- vn 
uA ^ vA 

Table 1: The opacity set included in the present simulation. Note that in neutral current 
reactions, represents all species of neutrinos {y^, i^ei^x) with Ux representing heavy- lepton neutrinos 
(i.e. Vn^fr and their anti-particles). 
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Striving for simplification of our modeling, we include a reduced set of neutrino-matter in- 
teractions (e.g., Table 1). Regarding the charged-current interactions, emission and absorption of 
electron neutrinos by neutrons (the first column in Tabic 1), emission and absorption of electron 
anti- neutrinos by proton (the second column), and emission and absorption of electron neutrinos 
by heavy nuclei (the third column), are included. For the neutral-current interactions, elastic scat- 
tering of all neutrino flavors off nucleons (the fourth and flfth column in Table 1), and the coherent 
elastic scattering (the sixth column) are included. For the cross sections of these reactions, we 
employ the ones summarized in Burrows et al. (2006) while omitting higher-order terms such as 
ion-ion correlations and weak magnetism for simplicity. 

In computing the neutrino cooling rate {Q^intr)^ furthermore consider the contribution from 
pair neutrino annihilation (5e-e+-i-j^p (Cooperstein et al. 1986), plasmon decay Qj^up (Ruffert et 
al. 1996), and nucleon-nucleon bremsstrahlung Qnn^nnvv (Burrows et al. 2006), which are also 
summarized in Sekiguchi (2010). Hence Q'^^intr can be expressed as, 

Qu,intr = Qi-+Qe-+Qi++Qe+ 

+ Qj^up + Qnn-^nnvp), (35) 

where and Q^-/+ represents the cooling rate by electron/positron capture on free nucleons 

and on heavy nuclei, respectively. 

Concerning the neutrino heating {Q^'^), wc only include the dominant heating reactions in the 
gain region, which is absorption of electron/anti-electron neutrinos by free nucleons. Then Q^'^ 
reads 

where co denotes neutrino energy, Ki^^i is the energy-dependent opacity for electron or anti-electron 
neutrinos (i.e. i = f g or i/g see Appendix A), and Juj-,T-Luj is the energy-dependent Eddington mo- 
ments, respectively. Yielding to the prescription of the so-called light-bulb scheme (e.g., Nordhaus 
et al. (2010)), a term of e"^"^ is introduced to vanish the neutrino heating smoothly as the opacity 
becomes higher inward down to the neutrino sphere. To take a gray approximation, we replace the 
energy integration in Equation (36) with the root-mean-squared (RMS) energy of the streaming 
neutrinos (e^^^j, see Appendix A for the deflnition) as 

j duK^,i{-J^X - K,i) - ^s.,i) {es.,ifh{-Ju>' - nn, (37) 

where h denotes the monochromatic opacity, in which the energy-dependence is replaced with the 
one of the rms energy (namely, Ki = hi ■ ■ Since J','H in Equation (37) is related to the 



'^The delta function in the left-hand-side of Equation (37) may be replaced by the Fermi-Dirac function. In the 
case, an additional factor of F2{'ri^, 0) (for the degeneracy limit -^2(0, 0) « 2) can be multiplied, which could potentially 
enhance the impacts of neutrino heating. 
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variables in the laboratory frame by Equations (28,29), the two-moment equations of E^, and Fi, 

with the source terms are finally closed. Having given explicit forms of Q^^, we are now moving on 
to summarize the GR hydrodynamic equations including the source terms in the next section. 



2.4. Hydrodynamic Equations 

From Equation (2), the hydrodynamic equations are written in a conservative form as, 

dtp, + di{p,v') = 0, (38) 
dtSi + djiSiV^ + ae^'t'Pdi) = -Sodia + SkdiP'' + 2ae^^Sldi(t> 

-ae'^'f'iSjk - Pjjk)da^''/2 - e'"t>aQ''jii„ (39) 
dtT + di{Sov' + e^'t'P{v' + P') - p,v') = + ae^'l'{Sij - P^ij)Av - SiD'a 

+e^'^aQ^ni„ (40) 

dtip^Ye) + di{p,Yev') = p.Fe, (41) 

where X = e^'^X, = pWe^'^, = /u* , f = Sq — p,, and Ye is the electron fraction. 

From Equations (32, 33), the source terms appearing in the right-hand-side of Equations (39, 
40) can be explicitly written as. 



HI 

Ui 



+ J2 e-^^-^{ss^fk,i-WF,i + Pl;,Uk) (42) 



w 



+ ^ e-^"-"{e,S~^A'^E, - fW)- (43) 

Fg in Equation (41) denotes the change in due to neutrino-matter interactions, which can be 
estimated in the same way as . Given an appropriate EOS, the hydrodynamic equations (38)-(40) 
are closed. We employ a tabulated EOS by Shen et al. (98) for heavy nuclei and uniform nuclear 
matter. Since Shen EOS contains contributions only from baryons, wc need to add contributions 
from electron/positron, and photon (see Appendix B for more details). 
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3. Initial Models and Numerical Methods 
3.1. Initial Models 

To assess GR and 3D effects on tlie neutrino-heating mechanism, we compute four models 
with a combination of SR or GR hydrodynamics in ID or 3D, which we label as ID-SR, ID-GR, 
3D-SR, and 3D-GR, respectively. We employ a widely used progenitor of a 15Mq star (Woosley 
&; Weaver (1995), model "sl5s7b2"). In our SR models, the space-time metric is assumed to be 
fiat (i.e., a = 1, = 0, jij = 6ij, (f) = 0, K^j = 0) and the self-gravity is solved by the Poisson 
equation, 

V^^NT = 47rS'o. (44) 

See Kuroda &; Umeda (2010) for more details, such as about the methods to treat self-gravity in 
the AMR structure, and also to evolve the space-time metric with GR hydrodynamics. 

To construct ID models in our Cartesian code, the following condition for the spacial velocity 
Ui (and also for Si) is enforced at every (hydro-) timestep, 

Ui = -riTx\ (45) 

As can be read, this operation eliminates the non-radial components of the flow velocity and 
momentum. Although the artificial elimination could potentially lead to the shift of the kinetic 
energy into the thermal one, our ID results (without and with neutrino heating/cooling) are in 
good agreement with the previous ID results as will be mentioned in Appendix C and section 4. 
This suggests that the manipulation is not severely bad for the sake of this study. 

The 3D computational domain consists of a cube of 10000^ km^ volume in the Cartesian 
coordinates. In our 3D models, we set the maximum refinement AMR level (Lamr, i-c, refine 
the AMR boxes in the vicinity of the center) at 5 in the beginning and then increment it as the 
collapse proceeds. We define the criterion to increment Lamr every time the central density exceeds 
1012,13,13.5 g cm"^ (see Kuroda & Umeda (2010) for more details). Each AMR level consists of 
8^ AMR blocks with a nested structure and each AMR block has 8^ cubic cells. It thus yields an 
effective resolution of Ax ^ 600 m at bounce in the center of our 3D models. 



3.2. Numerical Methods 

Since the hydrodynamic and transport equations (Equations (38)-(40) and (21, 22)) are all 
expressed in a hyperbolic form, they can be evolved by a standard high-resolution-shock-capturing 
scheme. We utilize the HLL (Harten-Lax-van Leer) scheme (Harten et al. 1983) to evaluate the 
numerical fluxes. A reconstruction of the primitive variables defined at immediate left /right of the 
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cell surface is performed by a monotonized central method (Van Leer 1977). The fastest (or right- 
going) and slowest (or left-going) characteristic wave speeds of fluid system, Xfiu, for i(G x,y,z) 
direction are obtained by solving the following second order equation, 

AX}i^ + 2B\fiu+C = (46) 

where 



A = 




B = 








C = 





f)2 + ^^^-aV\ (47) 



and Cs is the sound velocity (see Appendix B.3). 



Meanwhile, the fastest and slowest characteristic wave speeds of radiation system, Xrad, are 
assumed to have the same expression of the radiation pressure (Equation (23)) as 

-^rad ~ 2 ^rad,thin H 2 -^radjthick) 

(48) 

where Arad,thin and Arad,thick is determined by -P^^in -Pthick) respectively. According to Shibata 
et al. (2011), the fastest (slowest) wave speed in the optically thick or thin limit is evaluated by 
taking maximum (minimum) values, that is. 



' . 2TyV±\/a^y*(2VF2 + l) -2M^V ■ A 



for the optically thick limit (where = 7'^Uj/u*) and 



-/3' ±a^L=,-l3' + aE-^], (50) 



for the optically thin limit, respectively. With these wave velocities regarding the fluid and radiation 
component (Xfiu/rad), we deflne the HLL flux (Anton et al. 2006) as 

p A+F^ - A-Fr + \-\+{Qr - Ql) 

bHLL = = = , (51) 

A-f- — A_ 

where \ = X/a, L/R denotes the left/right states for the Riemann problem with F^/^ and Ql/r 
representing the advection and conservative terms, respectively. 

To ensure conservative laws at the interface of different AMR levels, we furthermore need to 
perform a "re/Zuxm^" procedure in estimating the numerical flux (see Kuroda & Umeda 2010). To 
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evolve the BSSN terms, we adopt the 4th order finite differencing for the spatial derivatives and the 
4th order upwind differencing for the advection terms (Zlochower et al. 2005; Etienne et al. 2008). 
Numerical tests are presented in Appendix C, in which we first show a ID adiabatic core-collapse 
test to validate the implementation of Shen EOS in the code, followed by the corresponding ID 
tests including neutrinos. 

4. Results 

First let us compare prebounce features among the four models (ID-SR, ID-GR, 3D-SR, and 
3D-GR) in section 4.1, and move on to the postbounce phase in section 4.2. Then in section 4.3, 
we will discuss the 3D/GR effects on the neutrino-heating mechanism. 

4.1. Infall and Bounce 

We begin our comparisons with the infall, bounce, and immediate postbounce phase. As seen 
from Figure 2, collapse to bounce takes slightly less time in our GR models (137 ms) compared to 
the SR models (141 ms), and the central density pc at bounce is approximately 2 times larger in 
the GR simulations than in the corresponding SR simulations (see the inset in Figure 2). For our 
non-rotating progenitor, the dynamics of collapsing iron core proceeds totally spherically till the 
stall of the bounce shock. This is the reason that the multi-D effects are invisible in the immediate 
postbounce phase. Hence we focus on the comparisons between the ID-SR and ID-GR model in 
the rest of this subsection. 

Figure 3 shows several snapshots of the lepton fraction (Itotai); electron fraction (Ye)-, and 
electron- type neutrino fraction (Yy^) for the ID-SR (left panel) and ID-GR model (right panel), 
respectively. After neutrino trapping (i.e. at a central density of a few 10^^ g cm^'^), the central 
lepton fraction (black lines) is shown to be conserved later on. In the trapped regions, the radial 
profile of the neutrino fraction (blue lines) is almost flat, while Y^,^ shows a gradual increase to 
satisfy the /3— equilibrium. 

From Figure 4, it can be seen that the lepton fraction at bounce (right end-point in density) 
is slightly larger for the ID-GR model (~ 0.364, red line) compared to the ID-SR model (~ 0.359, 
black line). The slight suppression of deleptonization is possibly because the neutrino opacity is 

effectively enhanced because of the more compact core in the GR model. Note that this trend 
is qualitatively in accord with the previous ID results in which a spectral neutrino transport was 
solved (e.g., Lentz et al. (2011); Miiller et al. (2010)). 

Figure 5 compares the mean energy of emergent neutrinos (e,y) between the ID-SR (left) and 
ID-GR model (right) (see Appendix A for definition). The mean energy is shown to be 20% larger 
(maximally near at bounce) for the GR model compared to the SR counterpart. This is also 
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because of the more compact core due to the GR hydrodynamics (Figure 6), leading to a more 
hotter neutrino sphere at smaller radii. 
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4.2. 3D hydrodynamics in the postbounce phase 

In the first 10 ms after bounce, the bounce shock turns into the so-called "passive" shock, 
which expands gradually outward with no positive radial velocities (Buras et al. 2006a) . As shown 
in Figure 7, the average shock radii until the passive expansion starts (tpb < 10 ms), remain almost 

close in all the models. It then diverges, which is more remarkable between the 3D models (red 
and blue lines). As will be discussed in the following, this is because multi-D effects (convection 
and SASI) sensitively affect the postshock hydrodynamic behaviors, also under the influence of the 
different treatment in gravity (SR versus GR). 

First of all, let us compare the shock evolution among the ID models in SR vs. GR (ID-SR 
(green line) or ID-GR (black line)). As expected, the shock radius is generally more compact for 
the ID-GR model (black line). On the other hand, the maximum shock extent is observed to be 
almost the same with each other ((-Rghock) ~ 130 km). Though we cannot unambiguously specify 
the reason, this trend was also seen in Miiller et al. (2010), who compared the shock radii in ID 
simulations with detailed neutrino transport between their GR model using the CFC approximation 
and the corresponding Newtonian model with the effective potential approach (e.g., right panel of 
their Figure 5). The maximum of (-Rshock) in Figure 7 for our ID models indicates the epoch when 
the passive expansion stops. Afterwards (tpb > 70 ms), the shock begins to shrink and a much more 
rapid recession is visible for the GR model (black line). The maximum shock radii and the shock 
recession timescalc obtained here arc again similar to those obtained in the previous ID results for 
the same progenitor model (e.g., Lentz et al. (2011); Miiller et al. (2010)). 

Four snapshots in Figure 8 are helpful to characterize the postbounce features in our 3D-GR 
model. The top left panel is for fpb ~ 10 ms, when the bounce shock stalls at a radius of ~ 90 
km (seen as a central blueish sphere). From the sidewall panels, the dominance of the £ = 4 and 
m = 4 mode can be seen in the postshock region, which is a numerical artifact inherent to the 
use of Cartesian coordinates. Comparing the top left to top right panel in Figure 8, the size of 
the outer sphere that marks the position of the shock (seen as greenish in the top right panel) 
becomes bigger, which is due to the passive expansion. At this stage, there forms the gain region in 
which neutrino heating dominates over neutrino cooling (e.g., Janka (2001)). The neutrino-driven 
convection gradually develops later on. The sidewall panels of the top right panel also indicate the 
growth of the postshock convection triggered by Rayleigh- Taylor instabilities. The entropy behind 
the standing shock becomes higher with time due to neutrino-heating, which can be inferred from 
a yellowish bubble in the bottom left panel. The high entropy bubbles {s[k b /haryon] > 10) rise 
and sink behind the standing shock. The shock deformation is dominated by unipolar and bipolar 
modes, which is a characteristic feature of the SASI. The size of the neutrino-heated regions grows 
bigger with time in a non-axisymmetric way, which is indicated by bubbly structures with increasing 
entropy (indicated by reddish regions in the bottom right panel). 

During our simulation time (100 ms after bounce), the shock radii can reach most further out 
for our 3D-GR model (red line in Figure 7). In contrast, the shock has already shown a trend 
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of recession in other models. Before we focus on the reason of it in the final section, let us next 
compare the activities of convection and SASI that we only touched on above. 

Figure 9 displays space-time diagrams of laterally averaged; (a) Brunt- Vaisala (BV) frequency 

wbv (top left panel) ; (b) the anisotropic velocity V]iniso(top right); (c) the pressure perturbation 
Ap in a logarithmic scale (bottom left), and (d) the net heating rate per baryon Qnct (bottom right) 
for our 3D-SR (top four panels) and 3D-GR models (the other four panels), respectively. Each of 
the quantities is defined as. 



WBV = sign(CL)Vl5eflfCL|, (52) 

where ^eflf represents effective gravitational potential that is estimated by taking a radial gradient 
of the potential in the Newtonian limit (^eff = d^Nx/dr). Ci, is the Ledoux criterion; 



dP 
ds 



ds dP 
— + 



(53) 



in which the neutrino contribution to entropy is taken into account where the /3-equilibrium is 
satisfied (Buras et al. 2006a). Following Takiwaki et al. (2011), Faniso is estimated as 

^aniso = \l{p[{Vr-{Vr)Y+vl+vl])/{p), (54) 

where {A) represents the angle average of quantity A. We define the normalized pressure pertur- 
bation Ap and the net heating rate per nucleon Qnet as 



Ap.ffi)!. (55) 



and 



Q„,.£!!^0!%, (56) 

P 

respectively. At first glance of Figure 9, one may not see any big differences between the 3D-SR 
(top four panels) and 3D GR models (the other), but indeed there are. Let us first discuss the 
properties of the four panels (a) - (d) taking the 3D-SR model as a reference and then proceed to 
focus on the differences between SR and GR. 

From panel (a) (top left) showing the BV frequency for the 3D-SR model, one can depict 
three typical convectively unstable regions in the postbounce phase; prompt convection (greenish 
region at tpb < 20 ms behind the shock^, postshock convection (seen as a narrow horizontal stripe 
behind the shrinking shock (just behind the outer most boundary labeled by white line), and PNS 



*Note that the shock is indicated by a white thin Une quickly rising after bounce and the passive shock stalls at 
a radius of i? ~ 150km . 
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convection (clearly seen as a thick horizontal stripe above the PNS at a radius of ~ 10 — 20 km 
later than ip^ > 60 ms). 

In our 3D results, the PNS convection develops only very weakly before ~ 60 ms postbounce. 
This is due to the stabilizing effect by a positive entropy gradient (see the positive gradient persisting 
outside the PNS surface (i? 10 km) in the right panel of Figure 10). Afterwards, the PNS 
convection gradually becomes vigorous with time as the negative lepton gradient nascent the PNS 
becomes remarkable (see the steepening slope of Yi near i? ~ 10 — 20 km in the left panel of Figure 
10). Comparing the black dotted line in Figure 10 (for ID-GR model at ipb = 60 ms) with the 
corresponding one (green line) for the 3D counterpart, the slope of the negative gradient is shown 
to become much smaller for the 3D model both in the profiles of lepton fraction (left panel) and 
entropy (right panel). This is a natural outcome of the convective overturns acting to wash out the 
local gradients. 

From panel (b) in Figure 10 showing the anisotropic velocity, the postshock convection (R > 
100km) is clearly seen as a reddish stripe running from top left to bottom right. Note in the 
panel that the prompt convection can be also seen like a narrow prolate spheroid colored by red 
f-t tph ^ 20 ms with 10 < i? < 60 km. As seen, convective overturns operate above the PNS 
(~ 10 — 20 km in radius) and below the shock ~ 100 km in radius. In-between, the region with 
smaller anisotropic velocity is formed (seen as a horizontal stripe colored by deep-blue at a radius 
of 30 — 50 km after tpb ~ 50 ms). By comparing to panel (d) (the net heating rate) to panel (c), 
the region is overlapped with the cooling region (Qnct < 0). The smaller anisotropic velocity there 
is because the infalling velocities in the cooling layer are so high that the convectively unstable 
material cannot stay there for long. Such a configuration has been already presented in 2D (Buras 
et al. 2006a) and 3D results (Takiwaki et al. 2011). 

Here let's see panel (c) not for the 3D-SR model but for the 3D-GR model for convenience. The 
accreting flows should receive an abrupt deceleration near at the bottom of the cooling layer (the 
dark colored region in panel (d) (bottom right)), below which the regions are convectively stable 
(panel (b)). There forms a strong pressure perturbation (seen as a greenish horizontal stripe near 
i? ~ 30 km in panel (c)). Subsequently the pressure perturbations are seen to propagates outward 
before they hit the shock (panel (c), seen as up-going stripes to the shock), probably leading to 
the formation of the next vortices. These features are akin to the so-called advectic-acoustic cycle 
(e.g., Foglizzo &; Tagger (2000); Foglizzo (2002); Scheck et al. (2008) and references therein), which 
is also observed in our 3D-SR model (top four panels). 

We now focus on the differences between the 3D-SR and 3D-GR models. Comparing the panel 
(a)'s between SR and GR, the unshocked core (regions below the PNS convection at tp^ > 50 ms) 
is shown to be more compact for the GR model. Between the pair models, Figure 11 compares the 
maximum of the pressure perturbation that the advecting vortices form near in the vicinity of the 
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deceleration regions^. As seen, the pressure perturbation after bounce is generally larger for the GR 
model (red line) compared to the SR model (black line) in our simulation time. This is presumably 
because stronger gravitation pull in GR makes the position of the coupling radius deeper, leading 
to produce more energetic acoustic waves. It is not straightforward to say something very solid only 
from the figure, but what we observed in our 3D-GR model (i.e. generation of stronger acoustic 
waves and the largest shock extent compared to the SR counterpart) does not seem, at least, 
unfavorable to drive neutrino-driven explosions. In the next section, we move on to discuss more 
in detail how 3D and GR would potentially impact on the neutrino-heating mechanism. 



Recalling that the neutrino heating rate can be symbolically expressed as oc L^{e1) (e.g., 
Janka (2001)), we first analyze the neutrino luminosities (L^) and the mean energies {{£u)) in the 
following. After that, we compare the dwell time to the neutrino-heating time in the gain region 
and discuss which one (3D-SR vs. 3D-GR) is most likely to satisfy the criterion to initial the 
neutrino-driven explosions. 

Figure 12 shows evolution of the neutrino luminosities of all the species (for v^, Vx (left panel), 
and i/g (right panel)) for all the computed models. Here the neutrino luminosity is calculated as 



where Q^' in Equation (32) takes into account all the cooling contributions. 

The spike in the luminosity corresponds to the so-called neutronization, when the shock 
propagates out through the v^. sphere. The peak v^. luminosity for the GR models is L^^ 3 x 10^^ 
erg s~^ (insensitve to ID or 3D), which is slightly luminous compared to those in the SR models 
{Ly^ ~ 2.9 X 10^'^ erg s"-*^). Using the same progenitor (Woosley & Weaver 1995), this trend is 
qualitatively similar to Bruenn et al. (2001). On the other hand, recent studies in which more 
detailed weak interactions are included in the Boltzmann transport have shown that the peak v^. 
luminosity becomes ~ 10% smaller for the GR models (e.g., Lentz et al. (2011); Miiller et al. (2010)). 
This may carry an important message that the Boltzmann transport should be implemented in the 
full GR simulations to obtain a ~ 10%-order accuracy, which is not small at all when speaking 
about the neutrino-driven mechanism. 

After the neutronization burst (tpb ^ 10 ms), the v^. luminosity for the GR models slightly 
increases later on, while it stays almost constant for the SR models during the simulation time 
(green and blue lines). The i>e luminosity after 50 ms postbounce (right panel in Figure 12) is 
highest for the 3D-GR model (red line), which is also the case for the Vx luminosity (left panel). 
Although the luminosities change with time, the luminosities generally yield to the following order. 



Scheck et al. (2008) termed it as the "coupling radius" in which the coupling of vortices and acoustic waves occur. 



4.3. 3D versus GR effects on the neutrino-heating mechanism 




(57) 
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for 3D-GR > IDGR, 3D-SR ~ ID-SR, 
for z^e, 3D-GR > IDGR, 3D-SR > ID-SR, 
for v^, 3D-GR > IDGR, 3D-SR > ID-SR. 

To summarize, both 3D and GR work to raise the neutrino luminosities in the early postbounce 
phase. As seen from the left panel in Figure 12, GR maximally increases the luminosity up to 
~ 50% (in 3D), while the maximum increase by 3D is less than ~ 20% (compare the Ve luminosity 
between the 3D-GR and ID-GR model). These results indicate that compared to the spacial 
dimensionality, GR holds the key importance to enhance the neutrino luminosities. 

Top two panels in Figure 13 compare the angle average of the RMS neutrino energy for (left 
panel) and (right panel) after the break-out burst (tpb > 10 ms). As seen, the RMS energies are 
highest for the ID-GR model (black line), followed in order by ID-SR, 3D-GR, and 3D-SR. In accord 
with the previous ID results (Lentz et al. 2011; Miiller et al. 2010; Liebendorfer et al. 2004; Bruenn 
et al. 2001), our 3D results (albeit limited to the early postbounce phase) support the expectation 
that the neutrino RMS energies increase when switching from SR to GR hydrodynamics. 

The reason for the higher neutrino energy is that the deeper gravitational well of GR produces 
more compact core structures, and thus hotter neutrino spheres at smaller radii. This is shown 
in the bottom panels in Figure 13 (compare the radii of the neutrino sphere between GR and 
SR models). The smaller neutrino energies for our 3D models compared to the corresponding ID 
models (top panels) is due to their larger neutrino spheres (bottom panels). In our 3D models, 
the shock expands much further out assisted by convection and SASI (e.g.. Figure 7), which also 
extends the positions of the neutrino spheres. The enlargement of the neutrino sphere in multi- 
D models is qualitatively consistent with the 2D post-Newtonian results by Buras et al. (2006a) 
including detailed neutrino transport. 

As mentioned above, GR increases the neutrino luminosities and energies, while the 3D hydro- 
dynamics works to make the neutrino energy smaller. What we like to discuss finally is whether the 
gain effects of GR could or could not overcome the possible loss effects of GR that should shorten 
the residency time of material in the gain region. And multi-D effects join in the game because 
they could potentially work against it to make the dwell time longer. 

A widely prevailing indicator to diagnose the onset of the neutrino-driven explosions is the 
ratio of the residency timcscalc ((tres)) to the neutrino-heating timescale ((theat)) in the gain region 
(e.g., Janka (2001); Thompson et al. (2005); Murphy & Burrows (2008))^°. To estimate (t^es), we 
employ the effective advection timescale (Equation (8) in Buras et al. (2006a)), in which (tres) is 
determined by the crossing time of mass shell between shock and gain radii. The local heating 



^'^If this indicator is greater than unity, i.e. (tres)/(theat) > 1, the neutrino heating proceeds fast enough to 
gravitationally unbind the fluid element, otherwise the matter is swallowed by the neutrino-cooling layer 
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timescale is estimated by the mass weighted average of the local heating timescale, 

^bind 

Theat = ~T j 
Wi^, total 

where we obey the Newtonian expression to estimate the local binding energy as 

and the net neutrino heating rate is calculated by Q^^ total = e^'^aQ^n^ (see, Equation (40)). In 
estimating the heating timescale, the numerical cells that satisfy both Ebind < and Q > are only 
taken into account. 

As seen from Figure 14, the shock revival seems most likely to occur for the 3D-GR model 
(red line) in our simulation time, which is followed in order by 3D-SR, ID-SR and ID-GR models. 
Thanks to a more degree of freedom, the residency timescale becomes much longer for the 3D 
models than for the ID models. In addition, the increase of the neutrino luminosity and RMS 
energies due to GR (Figure 13) enhances the timescale ratio up to the factor of ~ 2 for the 3D-GR 
model (red line) compared to the SR counterpart (blue line). Therefore our results suggest that the 
combination of 3D and GR hydrodynamics could provide the most favorable condition to trigger 
the neutrino-driven explosions. 

As expected from Figure 14, the shock revival will never occur afterwards for the ID models 
that have already shown the sign of a rapid shock recession. On the other hand, the curves for 
the 3D models stay constant for the last 30 ms before our simulation terminates. For the 15 Mq 
progenitor employed in this paper, the neutrino-driven explosions are expected to take place later 
than ~ 200 ms postbounce at the earliest (Bruenn et al. 2010) and it could be delayed after ~ 600 
ms postbounce (Marek &: Janka 2009) as already mentioned. The parametric explosion models have 
shown that the earlier shock revival is good for making the explosion energy larger (e.g., Nordhaus 
et al. (2010)). The onset timescale of the neutrino-driven explosions predicted in 2D models (Marek 
& Janka 2009; Bruenn et al. 2010; Suwa et al. 2010, 2011) could be shorter if the combination effects 
of GR and 3D would have been included. We anticipate that this can be a possible remedy to turn 
the relatively underpowered 2D explosions into the powerful ones. To draw a robust conclusion, 
the energy and angle dependence of the neutrino transport should be accurately incorporated in 
our full GR simulations with the use of more detailed set of weak interactions. This work is only 
the very first step towards the climax to investigate these fascinating issues. 
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5. Summary 

We presented the results from the first-generation multi-D core-collapse simulations in full GR 
that include an approximate treatment of neutrino transport. Using a Ml closure scheme with 
an analytic variable Eddington factor, we solved the energy-independent set of radiation energy 
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and momentum based on the Thome's momentum formahsm. To simplify the source terms of the 
transport equations, a methodology of multiflavour neutrino leakage scheme was partly employed. 
Our newly developed code was designed to evolve the Einstein field equation together with the 
GR radiation hydrodynamic equations in a self-consistent manner while satisfying the Hamiltonian 
and momentum constraints. An adaptive- mesh-refinement technique implemented in the three- 
dimensional (3D) code enabled us to follow the dynamics starting from the onset of gravitational 
core-collapse of a 15 Mq star, through bounce, up to about 100 ms postbounce in this study. By 
computing four models that differ by ID or 3D and by switching from SR to GR hydrodynamics, 
we studied how the spacial multi-dimensionality and GR would affect the dynamics in the early 
postbounce phase. Our 3D results support the anticipation in the previous ID results that the 
neutrino luminosity and the average neutrino energy of any neutrino flavor in the postbounce 
phase generally increase when switching from SR to GR hydrodynamics. This is because the 
deeper gravitational well of GR produces more compact core structures, and thus hotter neutrino 
spheres at smaller radii. By analyzing the residency to the neutrino-heating timcscalc in the gain 
regions, we pointed out that the criteria to initiate ncTitrino-drivcn explosions could be most easily 
satisfied in the 3D models, irrespective of the SR or GR hydrodynamics. Keeping caveats in mind 
the omission of energy- and angle-dependence of the radiation fields and the use of reduced set 
of weak interactions in the present algorithm, our results indicated that the combination of 3D 
hydrodynamics and GR should provide the most favorable condition to drive a robust neutrino- 
driven explosion. 

On top of the omission of the spectral and angle dependence of the neutrino transport, we 
think the most urgent task is to replace the leakage scheme with more realistic modeling of the 
source terms. In our 3D simulation, the numerical resolution behind the standing accretion shock 
is a few kilometers, which is not good enough to capture the growth of SASI accurately (Sato et 
al. 2009). The numerical viscosity is expected to be large especially in the vicinity of the shock, 
which may affect the growth of the SASI. It could also affect the growth of the turbulence in the 
postshock convectively active regions, which is very important to determine the success or failure 
of the neutrino-driven mechanism. To clearly see these effects of numerical viscosity, we need to 
conduct a convergence test in which a numerical gridding is changed in a parametric way (e.g. 
Hanke et al. (2011)). 

The most up-to-date neutrino transport code in core-collapse supernova simulations can treat 
the multi-energy and multi-angle transport in 2D (Ott et al. 2008) and even in 3D simulations 
(Sumiyoshi & Yamada 2012) but mostly in the Newtonian hydrodynamics (see, however, Miiller et 
al. (2011, 2012)). As was originally pointed out by Schwartz (1967) in the late 60's, our exploratory 
results also support the importance of GR to draw a robust conclusion to the supernova mechanism, 
indeed. The combined effects of GR and 3D^^ should affect not only the supernova dynamics, but 



^^It is worth mentioning that the MHD effects also remain to be studied (Kotake et al. (2004); Takiwaki ct al. 
(2004, 2009); Burrows et al. (2007a); Gullet et al. (2011); Kuroda & Umeda (2010); Obergaulinger & Janka (2011); 
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also the observational signatures of gravitational- waves (e.g., Miiller et al. (2011); Kotake et al. 
(2009a,b); Ott et al. (2011)), neutrino emission (e.g., Abbasi et al. (2011); Marek et al. (2009); 
Lund et al. (2010)), and explosive nucleosynthesis (e.g., Fujimoto et al. (2011); Thielemann et al. 
(2011)). Keeping our efforts to improve the caveats mentioned above, we are going to study these 
fascinating subjects one by one in the near future. 



A. Determination of Neutrino Spheres 

As illustrated in Figure 1, we have to calculate the neutrino optical-depth (tj^) to determine 
the position of the neutrino spheres. It can be done by solving the following differential equation 
(actually by a matrix inversion) 

= -i^v, (Al) 

r 

with an appropriate boundary condition (tu\j.^^ = 0). Here k^, represents the neutrino opacity 
of each species. Then the neutrino sphere is determined where the optical-depth exceeds 2/3 
(tj^ = 2/3). For the matrix solver, we take the same one to solve the Poisson equation (44). As 
shown in Table 1, we only include a reduced, but the most fundamental set of weak interactions 
in the supernova cores, which consists of the charged-current interactions; nve ^ &~Pi P^e ^ e~^n, 
UgA -H- e~A' and scattering processes; up o up, vn -H- vn, vA -H- vA. Note v in the scattering 
processes, represents all species of neutrinos (fg, Pg, v^). The opacity for electron, anti-electron, and 
heavy-lepton neutrinos can be expressed as 

Ki/e = l^ai'^e'n) + Ka{Vf,A) + Ks{Uen) + Ks{UeP) + Ks{Ue^), (A2) 



«Pe = l^ail^eP) + KsiUen) + Ks{Uep) + Ks{UeA), (A3) 

and 

«i/e = i^siv^n) + Ks{vxP) + tis{vxA), (A4) 

respectively (e.g., Ruffert et al. (1996)). Here the subindex a and s denote absorption and scattering 
processes, respectively. Detailed descriptions for the expressions of each opacity can be found in 
Bruenn (1985); Burrows et al. (2006). 

Since we do not transfer the number density of neutrinos in the present scheme, we evaluate the 
emergent root-mean-squared (rms) energy of neutrinos Sg^^i in the following way. Note here that 
the subscript i = Ue, denotes the neutrino species. We first project the neutrino sphere defined 
in the cartesian grids to the spherical polar grids, which gives us the position of the neutrino sphere 



Takiwaki & Kotake (2011), see also Kotake et al. (2006) for collective references therein). 
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expressed in the polar grids as Ru,i{d, 4>)- Then we identify Ss^^i with the energy of neutrinos at the 

neutrino sphere assuming that they stream freely outwards with possessing the information of the 
last scattering surface. Then Sg^ at arbitrary point {R, 9, 0) may be expressed as 

es^,i{R, 9, <j)) = e,,i{R,,i{9, 0), 6, </.). (A5) 

Here eu,i in the right-hand-side denotes the neutrino energy at Riy{6, (p), which is estimated by 

where is the Fermi-Dirac integral and iji, = /i^/kBT is the degeneracy parameter with /Xj/, T, 
ks representing the neutrino chemical potential, matter temperature, and Boltzmann constant, 
respectively. 



A.l. Neutrino diffusion terms 

Here we briefly summarize how to determine the diffusion term; 

Qy-^ffM^ey (A7) 

according to Ruffcrt ct al. (1996); Rosswog k. Liebendorfer (2003); Sekiguchi (2010). In the above 
equation, depending on k (=1 or 0), one can obtain the energy (or number) diffusion rate, in 
Equation (A7) and in Equations (32)- (33) are model parameters that affect the neutrino diffusion 
timescale and the position of the neutrino spheres, respectively. We adjust these parameters {an = 2 
and = 3/2) in such a way to fit the neutrino luminosity obtained in the ID results using the IDS A 
scheme (see Appendix C.2). {g^^ = gp^ = 1 and ^j,^ = 4) in Equation (A7) simply represents a 
multiplicity of neutrino species, rii, is the number density of neutrinos per each energy bin e,^ in 
thermal equilibrium with matter and expressed by the Fermi-Dirac distribution function as 



here /x,/ is the chemical potential of neutrinos. We define the diffusion time scale of neutrinos as 

Tdiff ^ 3^^r(£.) (A9) 



where Ax{£i,) is assumed to be Ax{ei,) = t{£v)/ n{£v)- We assume the optical depth and opacity 
can be expressed as 
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by neglecting the higher order correction terms of £1,. and are energy independent optical 
depth and opacity, respectively. Finally, the energy integration in Eq.(A7) is rewritten as 



_ 1 47rc f et ^ 



1 /] 'JT'f' 



where 771/ = fXi^/ksT is the degeneracy parameter of neutrino. 



B. Implementation of EOS 

We employ Shen EOS (Shen et al. 1998) based on the Thomas-Fermi approximation and a 
minimization of the free energy within a rclativistic mean field theory. The available data^^ is 
tabulated as a function of the three thermodynamic variables of density, temperature, and electron 
fraction as {p,T,Y(,). We smoothly interpolate/extrapolate the original data as; 10^'^ < p < 
10^^ g cm-^ with 200 equidistant intervals in a logarithmic scale, 10^ < T < 10^^ K with 200 
equidistant intervals in a logarithmic scale, and 0.01 < 1^ < 0.55 with 50 equidistant intervals in 
a linear scale. The interpolation is performed first by a bicubic interpolation in the p-T plane and 
then by a cubic interpolation for Yg direction. 

Thermodynamic variables such as the total pressure, internal energy, and entropy reads, 

P{p,T,Ye) = Pb + P^_+P^++P^, (Bl) 

e{p,T,Ye) = £6 + £e- +£e+ +£7, (B2) 
S{p, T, Ye) = Sb + Se- + Se+ + S7, (B3) 

where subscripts 6, e~, and e+ denote the contributions from baryon, electron, and positron, 
respectively. 



B.l. Supplement to the original Shen EOS 

Since the original Shen EOS contains contributions only from baryons, we need to add the 
remaining contributions from leptons and photons. Although it can be done straightforwardly by 
using formulae give in (e.g., Blinnikov et al. (1996)), we summarize them, for convenience, shortly 
in the following. 

To construct the leptonic EOS, we have only to determine the electron chemical potential p^- 
from a given data-set of proton fraction, density, and temperature (l^,/9, T). This can be done by 



'We use the updated version which is obtained from http://user.numazu-ct.ac.jp/ sumi/eos/#shen2011 
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the charge neutrality condition Yp = Ye, where Yg is defined by 



Y, 



rib 



(B4) 



where n^-/+ and = p/ruu is the number density of electrons/positrons and bayrons with 
being the atomic mass unit. n^-/+ can be expressed by 



«e-/+ 



(B5) 



where /3 = kBT/nieC? and r]~/'^ = fi^-/+/kBT with ks, me, T, /ig-/+ representing the Boltzmann 
constant, electron rest mass, temperature, and chemical potential of electron/positron, respectively 
(Blinnikov et al. 1996). -Ffc(r/, P) is the Fermi-Dirac integral, 



Fkiv, /3) 



oo 



+ l3x/2) 
e^-'J + 1 



1/2 



-dx. 



(B6) 



where the useful analytical formulae of the integral (with their derivatives) are given in Tooper 
(1969); Miralles &: van Riper (1996). Remembering that r/~ = — r/"'" is satisfied in the supernova 
cores due to high temperature (> 10^ K), one can find the solution of ;Ug- by Equations (B4) 
and (B5). Thus the total pressure, specific internal energy, entropy (per nucleon) can be readily 
calculated as 



3 H^tt'^ 



4^5 



/3'/'k/2(%,/3) + §i^5/2(%,/3) 



(B7) 



2 S 

-mtcr 



£e = V2-^/3-V2 F3/2(77e,/3)+^F5/2(r?e,/3) 



(B8) 



pee + Pe- n-ePe 

pTNAkB 

where A^^ is the Avogadro constant. 

The contribution from photons is expressed as, 

P =-oT^ 

1 7 ^ ' 

p 

~ 3pTNAkB ' 

where = 8Tr^k^/{15c^h^) denotes the radiation constant. 



(B9) 



(BIO) 
(Bll) 

(B12) 



-28- 



B.2. Primitive recovery 

Since we evolve hydrodynamic equations in a conservative form, we need to obtain primitive 
variables from the conservative ones. For the primitive recovery, we first solve the following si- 
multaneous equations to obtain Z = phW'^ and the Lorentz factor W (Cerda-Duran et al. 2008; 
Kuroda & Umeda 2010) 

{Z^ - S'^)W^ -Z'^ = (B13) 
T + D-Z + P{Z, W, YejYi) = (B14) 

for a given conservative set of variables (/9*, S^jt) and the electron/total lepton fraction Y^jYi ■ In 
the above equations, S"^ = ^"^^ SiSj and D = p^/e^'^ . P is the pressure and can be determined once 
the enthalpy h = Z/DW and the rest mass density p = D/W are given. We iteratively solve these 
equations by the Newton-Raphson method until the sufficient convergence is achieved. 

B.3. The sound velocity 

As given in Shibata Sz Sekiguchi (2005b). the sound velocity is expressed as, 



1 


dP 


P dP 
p^ de 

e 




Ti 


dp 


P_ 



(B15) 



where P and e include the sum of contributions from baryon, electron, and photon. Regarding the 
partial derivatives of the thermodynamical variables, we take a finite differencing of Shen's EOS 
table for the baryonic part, meanwhile we use analytical formulae of the Fermi-Dirac integrals given 
in Miralles & van Riper (1996) for the leptonic sector. 



C. Numerical Tests 

C.l. Core-Collapse tests with Shen EOS 

We first present the ID (in the same manner as ID-SR/GR models by neglecting the non-radial 
matter velocity and momentum) core-collapse run without neutrinos to validate the implementation 
of Shen EOS instead of the phenomenological one taken in the original code (Kuroda & Umeda 
2010). In the case of the adiabatic collapse, the so-called prompt explosion is expected to occur for 
the ISMq star (Woosley & Weaver 1995) as reported by Sumiyoshi et al. (2004) in their ID GR 
Lagrangian simulations using the same EOS. 

Figure 15 shows the profiles of density (left panel) and radial velocity (right panel) between 
the GR (solid line) and SR (dashed line) model, respectively. 
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As can be seen, the central density pc (left panel) and the infall velocity (right panel) becomes 
higher for the GR model, which bounces at = 4.5 x 10^^ g cm~^ with its inner-core baryon 
mass (Mic = O.QIM©) being ~ O.IM© smaller compared to the SR counterpart. After bounce, 
the prompt shock propagates through the entire iron core (left panel in Figure 16) for both of the 
models. For the GR model, the shock reaches at a radius of 1000 km at ~ 20 ms after bounce, with 
its explosion energy in the range of 1 - 1.5 xlO^^ erg (right panel in Figure 16), which is consistent 
with those obtained in Sumiyoshi et al. (2004). 

As for another test to check numerical accuracy, we show the violation of the average Hamil- 
tonian constraint Chm and the Arnowitt-Deser-Misner mass (ADM mass) Madm in figure 17. We 
adopt the following form for Chm (Shibata 2003) 



1 



hm 



Mba 



p^Ti 



+ 



e't'R 



+ 



27r(5o + E)e-'P\ + |^ (i^.i^^' - fif^) 



(CI) 



where Mbar is the proper rest mass and is the Hamiltonian constraint. 



=5</< 



n = D'Die'^ - ^ + MSo + E)e-'P + — AijA'^ - -K 



0. 



-^ADM can be written as 



(So + E)e-^ + — ( - - f^R,^e-^<t> 



(C2) 



(C3) 



The violation of the ADM mass (red line) here for our 3D-GR model becomes larger in the post- 
bounce phase, which is not surprising considering the complicated non-linear nature of the field 
equation and also the presence of the shock that makes the accuracy of the high order shock- 
capturing scheme down to the first-order scheme inevitably. However, Chm (black line) is generally 
kept less than 10^'^ in the postbounce phase, which seems to be slightly better than the typical 
ones obtained in the previous full GR simulations (e.g., Shibata (2003); Sekiguchi (2010)). 



C.2. Tests for Transport Scheme 

For numerical tests of our transport algorithm, we present the check for the trapped and 
streaming neutrinos, respectively. Note again that the sum of the streaming and trapped neutrinos 
is transported by the evolution equations (Equations (21), (22)). The streaming part can be 
estimated by subtracting the trapped contribution from the sum, because the trapped part can be 
simply determined by local hydrodynamic quantities (i.e. density, temperature, and Yg). 

Figure 18 shows comparison of the RMS neutrino energy of the trapped neutrino and the one 
obtained by the IDSA scheme below the neutrino sphere. As seen, the assumption of /3-equilibrium 
condition works well in the high density region and both of them show a quite similar profile there, 
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in which neutrinos are essentially trapped by matter. The agreement regarding the position of the 
neutrino sphere also certificates our estimate of the RMS neutrino energy (Equation A5). 

To check the properties of the streaming neutrinos, which is very relevant in the postbounce 
dynamics, wc check the following two points, which arc (1) whether the radiation energy flux falls 
with proportional to r^^ above the neutrino sphere and (2) whether the gain region, in which the net 
heating rate becomes positive, can be formed similar to previous studies. Figures 19 shows radial 
profiles of the radiation (energy) flux that is obtained by solving the closed set of the two-moments 
equations (Equations (21,22)). 

As seen, the energy fluxes change with r^^ (compare with the black line) irrespective of neutrino 
species outside the neutrino spheres. Note that the position of the neutrino spheres can be seen as 
the intersection point between the horizontal line (r = 2/3) and the solid lines. The radii of the 
neutrino spheres obeys a canonical order Ri^^ < Rp^ < R^^ . The emergent neutrino energy flux can 
be estimated by Airr'^, 47rr^ x {Fr,u^, Fr^p^, Fr^uJ ~ (9 x 10^^,8 x 10^^,4 x 10^^) erg s~^, which are 
all in good agreement with the luminosity defined by Equation (57) as plotted in Figure 12. 

Figure 20 shows evolution of the net heating rate and radial velocity along the x axis for our 
3D-GR model (see, Sec. 3) at selected postbounce epochs. As the passive shock propagates (from 
top left to bottom right panels), the gain region also gets larger. This reflects that the neutrino 
absorption on free nucleons predominantly takes place in the (enlarging) postshock region. The 
positive peak in the net heating rate is shown to be around 0.2GeV/nuc/s in the first 100 ms 
postbounce, which is in agreement with those in previous studies (e.g., Liebendorfer et al. (2001); 
Sumiyoshi et al. (2005); Ott et al. (2008)). 
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Fig. 1. — A schematic illustration how to model the source terms in our transport scheme. In 
this figure, trapped and streaming neutrinos are represented by 'Vtrap" and 'Vstream"; respectively. 
Qdiff Qintr^ ^^^^ Q^'^ denotes the coupling terms between them (see text for more details). 
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Fig. 2. — Central density pc as a function of time from initial collapse for our models of 3D-GR 
(red), ID-GR (black), 3D-SR (blue), and ID-SR (green), respectively. The inset is just zoom up 
near bounce, in which the time is measured from bounce (tph = 0) and the vertical lines represent 
the central density normalized by 10^^ g cm~^ in a linear scale. 




Fig. 3. — Profiles of the total lepton fraction (ItotaOi electron fraction (Y^), and electron-type 
neutrino fraction (Y^^) for the ID-SR (left panel) and ID-GR model (right panel) at times, when 
the central density reaches the value as indicated in the plots. 
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Fig. 4. — Comparison of the electron and lepton fraction versus central density during collapse for 
the ID-SR (black line) and ID-GR (red line) model, respectively. 
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Fig. 5. — Same as Figure 3, but for the profiles of the RMS neutrino energies for z^e (black) and 
(red). 
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Fig. 6. — Profiles of enclosed mass M{R) (solid) and the compactness parameter M{R)/R {dash- 
dotted) as a function of radius -R, in which lines are drawn every 2 ms in the first 10 ms postbounce 
for the ID-SR (red line) and ID-GR (black line) model, respectively. The ID-GR model has the 
maximum compactness parameter that is 10% larger, and the mass of the homologous core at 
bounce that is 20 % smaller compared to the ID-SR model. 




Fig. 7. — Evolutions of average shock radii as a function of post-bounce time tpb for the four variant 
models. 
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Fig. 8. — Three dimensional plots of entropy per baryon for four snapshots (top left; tpb = 10 ms, 
top right; tpb = 40 ms, bottom left; = 80 ms, and bottom right; tpb = 100 ms) for the 3D-GR 
model. The contours on the cross sections in the x = (back right), y = (back bottom), and 
z = (back left) planes are, respectively projected on the sidewalls of the graphs to visualize 3D 
structures. For each snapshot, the arbitrary chosen iso-entropy surface is shown, and the linear 
scale is indicated along the axis in unit of km. 
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Fig. 9. — Panels (a) to (d) display angle averaged; (a) Brunt- Vaisala frequency (wbv rns^^); (b) 
anisotropic velocity Vaniso normalized by 10^ cm s~^; (c) normalized pressure perturbation Ap (in 
a logarithmic scale) (d) net energy deposition rate per baryon Qnet [GeV nuc~^ s^^] for our 3D-SR 
(top four panels) and 3D-GR model (the rest four), respectively. Note that convectively unstable 
regions (i.e., wbv > 0) are only shown in panel (a) and the white line represents the contour of 
WBV = 0. 
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Fig. 10. — Profiles of angle-averaged total lepton fraction {left) and entropy {right) at times for 
the 3D-GR model, when the postbounce time is as indicated in the plots. Note that the result of 
ID-GR model at fpb = 60 ms (black dotted curve) is shown for comparison. 




Fig. 11. — Evolution of the maximum of the pressure perturbation Apmax in our 3D-SR and 3D-GR 
models. In taking the maximum, we set the radial range as 20 < i? < 50 km to cover the coupling 
radius (see panels (c) in Figure 9). 
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Fig. 12. — Neutrino luminosities of all neutrino flavors as a function of postbounce time (for v^, Vx 
(left panel), and for v^, (right panel), respectively. 




Fig. 13. — Evolution of the angle average RMS neutrino energy {upper panels) and the radii of the 
neutrino spheres (/ower panels) for (left) and v'e (right). Colors are as in Figure 12. 
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Fig. 14. — The ratio of the residency timescale to the heating timescale for the set of our models 
as functions of post-bounce time (see text for the definition of the timescales). 




Fig. 15. — Profiles of the rest-mass density {left) and the radial velocity (right) at times, when the 
central density reaches at 10^^'^^'^^ g cm~'^ (from the bottom up to the top in the left panel, each 
density is denoted by 12, 13, 14 and by "CB" at bounce in the right panel). Solid and dashed line 
is for the GR and SR model, respectively. In the right panel, the profiles of the sound velocity and 
the sonic point are indicated by black solid lines and black points at the intersection between the 
radial velocity and the sound velocity. 
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Fig. 16. — Evolution of the radial velocity (left panel) and the explosion energy (right panel) for the 
GR (solid line) and SR model (dashed line), respectively. "CB" (core bounce) and the postbounce 
time is shown for reference. Note that the explosion energy is defined in the Newonian limit that 
refers to the integral of the energy over all zones that have a positive sum of the specific internal, 
kinetic and gravitational energy. A smaller inner-core mass in the GR model leads to a smaller 
explosion energy because the amount of dissociation of iron nuclei becomes larger during the shock 
progagation. 




Fig. 17. — The violation of the Hamiltonian constraint Chm (black) and the ADM mass (red) in 
our 3D-GR model are plotted against the postbounce time. Refining procedures and enforcing the 
Hamiltonian constraint are done every time the number of the AMR blocks is increased by the 
AMR procedure (see section 3) and every time a restart of the simulation is needed such as after 
the maintenance of the supercomputing system. 
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Fig. 18. — Profiles of the RMS neutrino energy at times, when the central density reaches 10^^'^'^'^^ 
g cm~^ obtained by the present scheme (black solid line) and by the IDS A scheme (red triangle), 
respectively. 
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Fig. 19. — Radial components of neutrino energy flux F^^u in all flavors are plotted by color coded 
points {green Ue, blue i>e and red Ux). The profiles of optical depth Tj^ are shown by color-coded 
solid lines. The horizontal dash-dotted line is drawn for r = 2/3. As a reference, black solid line 
represents the slope of in the log-scale. The data are at 20ms after bounce in the 3D-GR model 
(see, section 3). 
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Fig. 20. — Profiles of the net heating rate Qnet[GeV by"-*^ s~^] {thick) and the radial velocity Vr 
(thin) along the x axis for our 3D-GR model at selected postbounce epochs. Note that Vr is 
normalized by 10^ km s~^. 



